Probability 2 (Data science Master, University of Lille) / CMF, (Centrale Lille, G3 SDIA)

---

Lab 1 - Discrete time homogeneous Markov chains

---

Guidelines (read carefully before starting)

Objectives: numerically simulate basic Markov chains (discrete time and discrete state space).

Setup: after retrieving the resources for the lab on moodle:

• place the .zip archive in a local folder (Computer -> Documents/Python/);
• unzip the archive .zip;
• rename the folder with the convention lab1_Name1_Name2
• duplicate the notebook file and rename it lab1_Name1_Name2.ipynb;
• [optional, possibly needed if working from Centrale's machines]
  • create a lab1 conda environment from the provided requirement.txt file

    conda create --name=lab1 --file=requirement.txt
    conda activate lab1
    # do not forget to deactivate the environment if needed
    # you can remove the environment once you are done
    conda env remove --name=lab1
    

  • launch jupyter notebook (the python environment will be the one from the conda environment lab1)
• at the end of the session, do not forget to transfer your work to your own network space if you are working on a machine from the school (do not leave your work on the C: drive).

Assessment → grade /20 (possibly converted later on to a grade ranging from F to A (A+))

This lab session will be evaluated, based on your answer to the exercises reported in a Jupyter notebook (e.g., this one) and any additional .py file produced. In particular:

• make sure the notebook you produce is properly annotated: each exercise should be introduced by a short sentence summarizing its context. Concisely answer each question from the guideline.
• relate the numerical results to the theory covered during the lecture whenever appropriate;
• codes without any explanations (in a text cell) will not be considered as a complete answer, and will only be attributed a fraction opf the grade attributed to the question.
• any code produced should be commented whenever appropriate;
• include appropriate axis labels and a relevant title to each figure reported in the notebook;
• document any function you introduce (using docstrings), and comment your code whenever appropriate (see, e.g., PEP 8 style recommendations).
  • use a reference docstring style, e.g., the google docstring style.
• give a clear structure to your notebook: longer theoretical explanations should be reported in markdown cells. Include LaTeX equations to improve the clarity of the document.

Additional evaluation criteria:

• Code clarity / conciseness / documentation
• Comment quality / overall presentation

Contents

• Exercise 1: Ehrenfest model
• Exercise 2: Simulation of a discrete time homogeneous Markov chain

import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import package as pack
import pandas as pd

---

Exercise 1: Ehrenfest model (↑)

Consider a system of $K = 30$ particles (labeled from 1 to $K$) evolving in a closed box. The box is divided into two compartments in contact with each other, respectively identified by an index, $0$ and $1$. A hole at the interface between the two compartments allows the particles to move from one compartment to the other.

Particle motion is modeled as follows: at each discrete time intant $n$, one particle is chosen uniformly at random and moved from its current compartment to the other. Let $X(n)$ denote the number of particles in compartment $0$ at time $n$.

1. Briefly justify that $\bigl(X(n) \bigr)_{n \in \mathbb{N}}$ is a Markov chain.

Pour décrire une chaine de Markov, on décrit son squelette et les transitions entre les différents états.

Avec $q = 0$

Ainsi, le modèle est une chaine de Markov.

2. Is the chain irreducible? (Positive) recurrent? Aperiodic or periodic? Prove each statement from your answer (recall the main steps covered during the exercise session).

La chaîne est irréductible car $\forall i,j$ $ \exists chemin $ tel que $P(i->j)>0 $

La chaine est récurente car on a pour $T_0$, temps de premier retour: $$ P(T_0<+\infty) = \sum_{n=1}^{\infty}P(T_0=n) = 1 $$

La chaine est périodique de période 2 car la chaîne peut être décomposé en de sous partie sur son univers : nombre impaire et nombre pair et c'est le plus petit nombre pour décomposer l'espace des états comme ceci. Si $X_n$ est pair (resp impair), alors $X_{n+1}$ est impaire (resp pair).

3. Recall the structure of the transition matrix, and encode it in Python (without any for loop).

Hint: use the function numpy.diag.

K = 30
p=1/2
diagonal = np.array([i/K for i in range(1,K+1)])
P = np.diag(diagonal,k=-1)+np.diag(1+1/K-diagonal,k=1)   # Définit la matrice de transition

4. Numerically verify that the binomial distribution $\mathcal{B} (K, 1/2)$ is invariant for the chain. This invariant distribution will be denoted by $\pi$ in the rest of the exercise.

Hint: you can use scipy.stats.binom.

B = binom(K, p)
x = range(K+1)
mu_bin = B.pmf(x)
plt.bar(x,mu_bin)
plt.title("Histogramme de la distribution binomiale")

Text(0.5, 1.0, 'Histogramme de la distribution binomiale')

On vérifie que µ est bien une mesure invariante avec la propriété : la mesure invariante est un vecteur propre à gauche pour P associé à la valeur propre 1

mu_bin.T@P - mu_bin

array([-5.16987883e-25,  3.30872245e-24,  7.94093388e-22, -3.81164826e-21,
        3.38813179e-20, -1.08420217e-19,  1.08420217e-19,  6.50521303e-19,
       -3.46944695e-18, -1.21430643e-17,  2.08166817e-17,  5.55111512e-17,
       -1.38777878e-16,  1.94289029e-16, -2.22044605e-16,  3.05311332e-16,
       -2.49800181e-16,  2.35922393e-16, -1.52655666e-16,  1.52655666e-16,
       -1.04083409e-16,  3.64291930e-17, -4.33680869e-18,  1.73472348e-18,
        0.00000000e+00, -2.71050543e-20,  3.72694497e-20, -1.27054942e-21,
        1.16467030e-21,  4.30133919e-23,  2.27474668e-24])

On vérifie ainsi de $\mu_{bin}$ est une mesure invariante.

5. Implement a Python function ehrenfest to simulate a trajectory of the chain for a system of $K$ particles, for initial distribution $\mu$ (for instance a Dirac centered in $0$, meaning that the compartment $0$ is empty at $n = 0$). The maximum number of time steps will be controlled by an input parameter $n_{\max}$. Random number generation will be controlled by a random number generator passed as an input to the function.

For an efficient implementation, do not use vector-matrix product with the transition matrix implemented in 3: rely on the description of the system instead.

Cette fonction est implémentée dans le fichier package.py

6. Simulate a trajectory of the chain starting in state $0$ for $n_{\max} = 5000$. Display the result in function of the time index $n$. Briefly describe the curve you obtained.

simulate_ehrenfest = pack.ehrenfest(mu=0,K=30,Generator = np.random.default_rng(),nmax = 5000)
plt.plot(simulate_ehrenfest)
plt.title("Trajectoire de la chaîne")
plt.xlabel("index de temps")
plt.ylabel("Indice de l'état courant")

Text(0, 0.5, "Indice de l'état courant")

7. Compare the empirical histogram of the trajectory obtained in 5. to the theoretical limit distribution $\pi$. What do you observe?

hist = plt.hist(simulate_ehrenfest,density=True,bins=K-5)
plt.xlabel("Etats visités au cours de la simulation")
plt.ylabel("Probabilité empirique")
plt.title("Histogramme des états visités au cours de la simulation des urnes d'Ehrenfest")

Text(0.5, 1.0, "Histogramme des états visités au cours de la simulation des urnes d'Ehrenfest")

On observe à la limite la distribution semble décrire la distribution de la loi binomiale.

plt.hist(simulate_ehrenfest,density=True,bins=K-5,color="red",label="Histogramme de densité des valeurs simulées")
plt.bar(x,mu_bin,color = "grey",label="Histogramme loi binomiale")
plt.legend()
plt.show()

8. a) Modify the function defined in 1. so that it returns the return time to state 0, defined as $T_{0,0} = \inf \bigl\{ n > 0, X(n) = 0 \mid X(0) = 0 \bigr\}$.

print(f"Le temps de retour pour un K=10 est de (simulation sur un trajet): {pack.return_time_0(K=10,Generator = np.random.default_rng())}")

Le temps de retour pour un K=10 est de (simulation sur un trajet): 110

8. b) [Optional] Run several chains (about 5, ideally in parallel) for $K = 10$, $n_{\max} = 5000$, and compare the empirical average of $T_{0,0}$ to $\pi(0)$. What do you observe?

Hint: a good tutorial showing how to run functions in parallel in Python is available here.

import multiprocessing as mp print("Number of processors: ", mp.cpu_count()) # Step 1: Init multiprocessing.Pool() pool = mp.Pool(mp.cpu_count()) # Step 2: `pool.apply` the `howmany_within_range()` results = [pool.apply(return_time_0, args=(10,np.random.default_rng())) for i in range(10)] # Step 3: Don't forget to close pool.close() print(results)

T_list = [pack.return_time_0(K=10) for i in range(2000)]
plt.hist(T_list,bins=100,density=True)
plt.title("Temps de retour en 0 pour K=10 avec 2000 simulations")
plt.xlabel("Temps de retour en 0")
plt.ylabel("Probabilité de ce temps de retour")
plt.show()

L_k=[]
for K in range(1,20):
    T_list = [pack.return_time_0(K) for i in range(2000)]
    L_k.append(np.mean(T_list))

plt.plot(range(1,20),L_k)
plt.xlabel("Evolution du paramètre K")
plt.ylabel("temps de retour en 0 moyen ")

Text(0, 0.5, 'temps de retour en 0 moyen ')

8. c) Comment on the possibility of numerically observing the chain returning to its initial state as $K$ increases.

...

---

Exercise 2: Simulation of a discrete time homogeneous Markov chain (↑)

Let $\bigl( X(n) \bigr)_{n \in \mathbb{N}}$ be a discrete time homogeneous Markov chain defined by the following initial distribution $\mu$ and transition matrix $P$

$$ \mu = [0, 1, 0, 0, 0, 0], % \quad % P = \begin{pmatrix} 1/2 & 1/2 &0 &0 &0 &0 \\ 1/4 &0 &0 &1/4 &1/4 &1/4 \\ 1/2 &0 &0 &0 &0 &1/2 \\ 0 &1/2 &0 &0 &1/2 &0 \\ 0 &1/3 &0 &1/3 &0 &1/3 \\ 0 &1/3 &1/3 & 0 &1/3 &0 \end{pmatrix}. $$

1. What can you say about the Markov chain $X$? (irreducibility, positive recurrence, periodicity, ...). Justify each of your claim, citing the relevant results from the lecture.

Ici les propriétés de la chaine de markov sont directes, nous détailleront peu les preuves :

• La chaine est irréductible car il existe une unique classe d'équivalence notée $C_{eq} = {0,1,2}$ En effet pour chaque couple i,j d'états, $\exists$ un chemin de $i->j$ avec $P(i->j)>0$ Ainsi la probabilité d'aller de i en j en un nombre d'état fini est non nulle.

• Cette classe d'équivalence est fermée donc elle est récurrente (une classe d'équivalence à soit tous ces états récurrents, soit tous ces états transitoires), donc tous les états sont récurrents.

• Cette chaine est de période 1 car il n'existe pas de séparation de l'espace d'état en sous espace tel que la probabilité de passer d'un sous-espace à un autre soit de 1. (ie : Que l'on soit dans l'état 0, 1 ou 2 il est possible d'aller directement dans les deux autres états.)

• La chaine étant irréductible, il existe une unique mesure invariante (donc proba invariante en normalisant). Ainsi la chaine est positive récurrente.

2. Write a function simulate_dthmc simulating the trajectory of the Markov chain $\bigl( X(n) \bigr)_{n \in \mathbb{N}}$ for $n_{\max}$ time steps. The signature of the function should include the following elements:

• list of inputs: transition matrix $P$, initial distribution $\mu$, number of time steps $n_{\max}$, random number generator;
• output: array of lenght $n_{\max}$ representing the trajectory.

To this end, do not use matrix vector products, which would lead to an extremely inefficient algorithm in this case.

P = np.array([
        [0.2,0.7,0.1],
        [0.9,0,0.1],
        [0.2,0.8,0]])
distribution_init = [0,1,0]

3. Simulate a trajectory of the chain for $n_{\max} = 2000$ starting from $X(0) = 1$. Plot the histogram of the states visited by the chain.

nMax=2000
states = pack.simulate_dthmc(transitionMatrix=P,distribution=distribution_init,nMax=nMax)
states=pd.DataFrame(states.astype(int),columns=['Etat'])

pack.show_histrogramm(states=states)

4. Determine numerically an invariant distribution $\boldsymbol{\pi}$ of the chain (e.g., based on an eigendecomposition numpy.linalg.eig). Is it unique? Compare it to the histogram obtained in 2 (superimpose graphs). What can you conclude?

La chaine étant récurrente et l'espace d'état fini, nécessairement la mesure invariante est unique. En effet d'après le Th de Perron-Fröbenius : Pour $P$ matrice stochastique et irréductible, 1 est valeur propre de cette matrice $\exists !$ vecteur à coefficients positifs associé à cette valeur propre tel que $\sum_{i=1}^{S}\mu_{i} = 1$

w,v=np.linalg.eig(P.T)

pack.show_superImposed_Histogramms(states,v)

Vecteur propre associé à la valeur propre 1 :           [-0.75529182 -0.64035611 -0.13956479]
Vecteur propre associé à la valeur propre 1 normalisé : [0.49197861 0.4171123  0.09090909]

Nous constatons en effet que la probabilité de visite de l'état 2 est autour de 0.09 contrairement à celle de visiter l'état 1 ou 0 qui vaut 0.41 ou 0.49 Ainsi en simulant la chaine un grand nombre de fois, nous devrions tomber sur cette probabilité invariante.

5. a) Compute $\mu_n \triangleq \mu P^n$, the probability distribution of $X(n)$. What is the limit of $\mu_n$ as $n$ goes to $+\infty$? Illustrate the result numerically.

nStep=50
distribution_n,history = pack.compute_distrib_n(P,distribution_init,nstep=nStep)

pack.show_histogramm_vector_df(pack.convert_vector_to_df(distribution_n))

$\mu_{n}$ tend vers la distribution invariante de la chaine de Markov, l'histrogramme ci-dessus illustre bien cette convergence vers la probabilité invariante de la chaine.

5. b) Display on the same graph the curves $n \mapsto \mu_n(i)$ for $i = 1, \dotsc , 6$, and compare with $\pi$. Display on another graph the function $n \mapsto \Vert \mu_n - \pi \Vert_1$, where $\Vert \cdot \Vert_1$ is the $\ell_1$ norm. What does each of these curves illustrate?

history = np.array(history)
history= history.reshape(50,3)
plt.figure()
plt.plot(history[:,0],label="Evolution de la distribution en 0")
plt.plot(history[:,1],label = "Evolution de la distribution en 1")
plt.plot(history[:,2],label = "Evolution de la distribution en 2")
plt.xlabel("Nombre de steps")
plt.legend()

<matplotlib.legend.Legend at 0x15bd7df7ac0>

pi_vect = np.array([pack.norm_eigen_vector(v) for i in range(0,nStep)])
diff = np.linalg.norm(history - pi_vect,ord=1,axis=1)
plt.figure()
plt.plot(diff)
plt.title("Evolution de l'écart entre la mesure invariante et son estimation")
plt.xlabel("Nombre de steps")
plt.ylabel("Norme l1 de la différence")
plt.ylim(0,1)
plt.show()

Ces résultats montrent que la distribition des probabilités de chaque état converge ici très rapidement. (1er graphe) Le second graphe montre l'écart entre cette distribution théorique et celle estimée.

6. For each state $i \in \{1, \dotsc, 5 \}$, simulate 100 trajectories starting from the state $i$ until the return time to $i$. For each state, compute the (empirical) average return time. Compare with its theoretical value.

MoyennePerState=[]
for state_considered in range(0,3):
     states = pack.simulate_dthmc(transitionMatrix=P,distribution=distribution_init,nMax=nMax)
     positionOfReturn = []
     for step  in range(len(states)):
          if states[step]==state_considered:
               positionOfReturn.append(step)
     MoyennePerState.append(pack.computeMoyenneReturn(positionOfReturn))

print(f"Le nombre moyen de retour pour chaque état est défini par (moyenne empirique)")
print(f"Moyenne du nombre de step pour le retour en 0 : {round(MoyennePerState[0],2)}")
print(f"Moyenne du nombre de step pour le retour en 1 : {round(MoyennePerState[1],2)}")
print(f"Moyenne du nombre de step pour le retour en 2 : {round(MoyennePerState[2],2)}")

Le nombre moyen de retour pour chaque état est défini par (moyenne empirique)
Moyenne du nombre de step pour le retour en 0 : 2.04
Moyenne du nombre de step pour le retour en 1 : 2.36
Moyenne du nombre de step pour le retour en 2 : 11.27

En théorie : Pour une chaine de markov positive récurrente irréductible et apériodique $lim_{n \to +\infty}P(X_{n}=i) = \mu _{i} = \frac{1}{E(T_{i})}$ Avec µ la mesure invariante pour le système

On en déduit les valeurs théoriques des temps de retour :

tabTempsRetour = [1/val for val in pack.norm_eigen_vector(v)]
print(f"Moyenne du temps de retour en 0 : {round(tabTempsRetour[0],2)}")
print(f"Moyenne du temps de retour en 1 : {round(tabTempsRetour[1],2)}")
print(f"Moyenne du temps de retour en 2 : {round(tabTempsRetour[2],2)}")

Moyenne du temps de retour en 0 : 2.03
Moyenne du temps de retour en 1 : 2.4
Moyenne du temps de retour en 2 : 11.0

Il y a évidemment un écart du à notre estimation en temps fini mais il est clair que notre estimation converge bien vers le nombre théorique de retours.

import os
os.system('jupyter nbconvert --to html lab1_notebook.ipynb')

0